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We generate equilibrium configurations for the three and four dimensional Ising spin glass with 
Gaussian distributed couplings at temperatures well below the transition temperature T c . These 
states are analyzed by a recently proposed method using clustering. The analysis reveals a hierar- 
chical state space structure. At each level of the hierarchy states are labeled by the orientations of 
a set of correlated macroscopic spin domains. Our picture of the low temperature phase of short 
range spin glasses is that of a State Hierarchy Induced by Correlated Spin domains (SHICS). The 
complexity of the low temperature phase is manifest in the fact that the composition of such a spin 
domain (i.e. its constituent spins), as well as its identifying label, are defined and determined by 
the "location" in the state hierarchy at which it appears. Mapping out the phase space structure by 
means of the orientations assumed by these domains enhances our ability to investigate the overlap 
distribution, which we find to be non-trivial. Evidence is also presented that these states may have 
a non-ultrametric structure. 

I. INTRODUCTION 

Whereas equilibrium properties of infinite range jjj spin glasses are completely understood within the framework 
of replica symmetry breaking (RSB) spin glasses with short range interactions are the subject of considerable 

current debate and controversy. Open questions address the nature of the low temperature phases and their 

theoretical description. Resolution of these issues by experiments or simulations is hindered by the extremely long 
relaxation time required for equilibration. 

The most widely studied model of a short-range spin glass is the Edwards- Anderson model of an Ising spin glass 

m 

where (ij) denotes nearest neighbor sites of a simple (hyper) cubic lattice in D dimensions (we will consider D — 3 
and D — 4) with periodic boundary conditions, Si = ±1, and the couplings, Jij, are independent random variables 
taken from a given distribution. The most commonly studied distribution, and the one we study here, is a Gaussian 
distribution with zero average and standard deviation J = 1. 

The high temperature phase of the model is a disordered paramagnet. As the temperature decreases below a critical 
temperature T c , the system (in three or more dimensions) undergoes a transition into a frozen spin-glass phase. In the 
spin glass phase, phase space is divided into "valleys" which we define as as an ergodic subset of the phase space, i.e. 
a maximal subspace that the system can span (or visit) as the time tends to infinity. For a finite system the definition 
is less clear, but a valley is usually referred to as a part of the phase space surrounded by free energy barriers, whose 
height diverges as the system size L — > oo. 

This definition of "valley" may not be identical to the notion of a "pure state" which has been used extensively in 
the literature [0-0 ■ and which is defined in terms of the set of correlation functions in a fixed finite region inside 
the system as L — > oo with some specified boundary conditions. In particular, it was recently emphasized [l^] 
that a spin glass can in principle have many thermodynamically important valleys but just two pure states. This is 
realized when there are many valleys with free energies which differ by an amount of order unity, and configurations 
taken from different valleys have a vanishing density of (relative) domain walls as L — > oo (a domain wall is a surface 
separating a region where the two configurations are identical from a region where they are opposite). In contrast, if 
the density of domain walls is finite {i.e. the domain walls are space-filling), there is a non- vanishing probability to 
have a domain wall in any finite region of the system, and thus to have more than two pure states. In this paper we 
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will be mainly concerned with the number and organization of valleys, and we will not investigate whether multiple 
valleys correspond to multiple pure states as defined above. In the following, by state we will always mean a microstate 
or spin configuration. 



A. RSB, droplet and TNT scenarios 



There are two traditional pictures of the spin glass phase; the droplet picture and RSB. According to the droplet 
picture of Fisher and Huse , the low energy excitations are in the form of droplets - compact regions with low 
surface tension that flip collectively. For a droplet of size L the typical (e.g. median) free energy Fl scales as L , 
where 9 is a dimension dependent exponent. Furthermore, the surface of these excitations has a vanishing density for 
large L. Therefore, thermodynamically important configurations have a vanishing density of relative domain walls, 
and hence a trivial overlap (defined below) over any finite region. It follows that in this approach within any finite 
region there are only two pure states, related by spin-flip symmetry. 

A parameter commonly used to measure domain wall density is the link overlap and its distribution. Denote a con- 
figuration (or state) of an TV-spin system by S M = (Sf , Slf,..., S%f). The link overlap q^" k between two configurations 
S^ 1 and S" is defined by 

^ = ^£WW' ( 2 ) 
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where the sum is over pairs of neighbor sites and 7-/V is the number of bonds in the system. If the domain wall density 
vanishes, then the distribution P(q hnk ) of the link overlap will be trivial: P(q hnk ) = 5(q hnk — qo). At T = one has 
qo = 1, while </o decreases for T > and becomes zero at T c . 

Another parameter commonly considered is the spin overlap q^ between configurations S M and S"; 
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If there are only two pure states, as in the droplet model, the local overlap distribution, obtained in a finite part of 
an infinite system, would be trivial for all T < T c , i.e. P(q) = 0.b[d(q — qEA) + 5(q + Qea)], where qEA is the average 
overlap inside a pure state. In addition, most conventional interpretations of the droplet picture Jl6|,[l7j argue that 
the global P(q), obtained from overlaps over the whole system, would also be trivial. This is realized if the droplets 
(with positive 0) are the only relevant excitations over all length scales. However the work of Huse and Fisher [Q, 
and also Newman and Stein Hjl8| , is formulated in a sufficiently general fashion to accommodate a non-trivial global 
P(q) if this arises from multiple valleys with non space-filling domain walls. In this situation, one would have a trivial 
link overlap distribution P(q hnk ) in the infinite system size limit. Even though the global P(q) would be non-trivial, 
the local P(q), would be trivial because a vanishing density of domain walls means that the probability that a domain 
wall goes through a fixed finite part of the infinite sample also vanishes. 



Numerical work has, so far, indicated a non-trivial global P(q) [[19 20 1. For example, Marinari et al. |2C]] have 
used parallel tempering (|^J2^] to sample 3D Ising spin glasses of sizes up to L = 16 and for temperatures down to 
T = 0.7 ~ 0.74T C . They have found that P(q) is non-trivial, and P(0) does not vanish. 

In the RSB picture, the Parisi theory, which is exact for the infinite range model JjJ, is assumed to also apply 
to short range systems. Within the RSB solution, both P(q) and P(q llnk ) are non-trivial for < T < T c , which implies 
that the system has many valleys and also many pure states. RSB suggests a tree-like hierarchical structure for the 
pure states. At every level of the hierarchy the states are divided into sets, so that the states in a given set are closer 
to each other than to states in other sets. At the next level down the hierarchy these sets are divided into subsets, 
and so on. Furthermore, according to the RSB solution the distances between the pure states exhibit ultrametricity 
||: the overlap between any two states is determined only by the lowest level in the hierarchy, at which they still 
belong to the same set. This means that for any triplet of pure states /i, v and p the following relation always holds: 

Qfiu > min^p, q vp ) . (4) 



Recently, a mixed picture has been proposed on the basis of numerical results of ground state computations [|12|-[14| , 
in which P(q) is non-trivial but P(q llnk ) is trivial (hence referred to as TNT; for Trivial and Non- Trivial). Houdayer, 
Krzakala and Martin [^JT^ I demonstrated the existence of macroscopic excitations with low energy cost in 3D Ising 
spin glasses of sizes up to L = 11. This suggests that the spin overlap distribution, P(q) is non-trivial at finite 
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temperature. Their results also indicate that the surface of these excitations is not space-filling, which suggests that 
the link overlap distribution, P(q hnk ), is trivial. 

Palassini and Young Q| studied changes to the ground state of a spin glass when a weak perturbation is applied 
to the bulk of the system. They considered short range models in three and four dimensions as well as the infinite 
range SK model and the Viana-Bray model. The results for the SK and Viana-Bray models agreed with the replica 
symmetry breaking picture as expected, but the data for the short range models agreed with the TNT picture. Effects 
of the type of perturbation considered in Ref . JL4| on RSB have been investigated by Franz and Parisi |23| . 

Katzgraber et al. measured directly the distributions P(q) and P(q hnk ) at finite temperature using parallel 
tempering |21 22 Monte Carlo, for 3D systems of linear size L < 8 at temperature T > 0.2, and 4D systems with 
L < 5 and the same temperature range. Extrapolating their results to large sizes they found that the variance of 
P(q lmk ) vanishes as L — > oo, and the distribution converges to delta function. They also found the distribution P(q) 
to be non-trivial, as in 19 2(J, so their results also agree with the TNT picture. In the TNT scenario there are many 
valleys separated by free energy barriers, but only two pure states p5|jL5 |. 

Although several pieces of work Jl2| [Ti 24 , 2^] supported a vanishing density of domain walls (and hence a fractal 
dimension of the domain walls, ds, less than the space dimension), a large extrapolation is involved in deducing this 
result, and Marinari and Parisi [E7W29] have argued, based on their own data and a somewhat different analysis, that 
actually d s — D, which corresponds to RSB. 



B. SHICS: State Hierarchy Induced by Correlated Spin Domains 



Very recently a new method of analysis of the structure of the low temperature phase of short range spin glasses has 
been introduced jj(J[n]. Evidence for a novel picture of this phase, which is consistent with the TNT scenario, but 
inconsistent with RSB (since there is no ultrametricity), has been presented J|(J on the basis of a "clustering analysis" 
of the degenerate ground states of the model (|l|) with = ±1 couplings. We denote this by "State Hierarchy Induced 
by Correlated Spin Domains" (SHICS). 

In this picture there is a hierarchical tree-like structure of the states as in the RSB solution. The highest levels 
of the state hierarchy, are schematically illustrated in Fig. [|. At the first level of hierarchy the states divide into 
sets C and C, such that a state in C has a counterpart with the same energy in C, obtained by flipping all the spins. 
This equality of the energies follows, of course, from the symmetry of the Hamiltonian in zero field. However, this 
symmetry information is not imposed on the analysis; the method finds it by itself. In fact, it is not trivial, for a 
spin glass, to divide the states into two clusters such that every state in C has its reversed state in C. Suppose, for 
example, that one has two states fi and v, and states \x and v with reversed spins, such that the spin overlap q^ v is 
close to zero. Should one put v or v in the same cluster as /i? The analysis, used in Ref. []30f and here determines 
which one it is. 

Many of the spins stay, with high probability, in the same relative orientation in most of the states C. Most of 
these form a contiguous cluster Q\, see Fig. [I]. Among the remaining spins, an apparently macroscopic fraction form 
a contiguous domain, Q2, such that the spins in it maintain, with high probability, their relative orientation in nearly 
all the states of C. Hence C divides into two sub-clusters of states, C\ and C2, depending on the orientation of Q 2 , 
see Fig. |l|. In general, the domains Q\ and Q2 are distinct. In many samples, further levels of the hierarchy, with 
successively smaller domains Q3, ■ ■ ■ can be clearly resolved, as discussed later. The excitations obtained by flipping 
the domains Ga, ■ • • appear to correspond to the large scale, low energy excitations investigated by Krzakala and 
Martin and Palassini and Young JlJ] . Note that the local (or link) overlap was not investigated in Ref. ||(| . 

By contrast, in the conventional interpretations of the droplet picture |16|l7j], the only substantial division of 
the states would be into C and C, and any further divisions emerging from the analysis would only correspond to 
microscopic spin domains. In the RSB scenario there would be a hierarchical structure to the states, similar to what 
we find here, but the nature of the spin domains would appear to be different, see e.g. [Q. We will discuss these 
differences further in Sees. IV and VII. 

The purpose of the present paper is to use use the methodology of Ref. Q to investigate whether the same picture 
of the spin glass phase found there also occurs for a spin glass with Gaussian couplings (which has a unique ground 
state apart from spin reversal) at finite temperatures. Both three and four dimensions are studied. We find that our 
data do fit this picture quite well. We also present here full details of the method. 

Readers who like to skip ahead will find the picture of state clusters and spin domains that were obtained at T = 0.2 
for a particular bond realization, conveniently summarized in Fig. [l2| The corresponding overlap distribution P(q) is 
presented in Fig. [l4| (a). 

The numerical procedure and parameters that were used in our simulations are described in Sec. ||. In Sec. |lll| we 
present the clustering methodology which we use to identify the states hierarchy, as described in Sec. |v| In Sec. |v| 
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FIG. 1. Schematic representation of the SHICS picture; the two largest spin domains and the first two levels in the hierarchical 
organization of the states are shown. The structure of the states is explained by the spin domains' orientations; e.g. in the 
states of the two sets Ci,Cz, the spins of the larger domain, Q\, have the same orientation, whereas the spins of the smaller 
domain, @2, have flipped. Spins not in Q\ or Q2 are in smaller domains which are not resolved at this level of the hierarchy. 



we use the hierarchical partition of the state space to obtain the spin domains, show that their sizes scale with the 
system size and their correlation does not approach unity as L — > 00. We also show that these spin domains, that 
were identified on physical grounds, can also be obtained by a cluster analysis of the ./V spins. Those domains yield a 
non-trivial overlap distribution P(q) with peaks corresponding to the different domain sizes, as we show in Sec. VI. 
Since we find that the average correlation between spins in different domains does not approach unity with increasing 
system size, P(q) will remain non-trivial as L — > 00. The nature of our pictu re appears to yield a non-ultrametric 
structure, as indicated at the end of Sec. |iv| and demonstrated in Sec. VI] , in which we present a para meter for 
ultrametricity and measure its distribution. Finally, our method and findings are summarized in Sec. VIII . 

After this work was completed we received a preprint from Marinari et al p3| who have adopted and adapted the 
methodology of pOpl| ] to study the = ±1 model in d = 3 at a single temperature (T = 0.5 - whereas here we 
considered T = 0.2 and T = 0.5 for the Gaussian model). They also confirmed that the previously observed SHICS 
scenario |3C[ ] of a tree-like structure of the states, governed by correlated spin domains, remains valid at a non-zero 
temperature. 



II. NUMERICAL METHOD 



We simulate the Hamiltonian in Eq. ([!]) using the parallel tempering Monte Carlo method pl] , p2| . In this tech- 
nique, one simulates several identical replicas of the system at different temperatures, and, in addition to the usual 
local moves, one performs global moves in which the temperatures of two replicas (with adjacent temperatures) are 
exchanged. This greatly speeds up equilibration at low temperatures. The detailed balance condition for temperature 
exchanges is satisfied by accepting these moves with probability min [exp(A£A/3), 1], where AE = — E v , E^ and 
E u are the (total) energies of replicas [i and z/, and A/3 = /3 M — f3 u is the difference in inverse temperatures. 

We choose a set of temperatures Ti,i = 1,2, ■■ • ,Nt, in order that the acceptance ratio for the global moves is 
satisfactory, typically greater than about 0.3. We use the test for equilibration discussed in Rcf. 24|, which involves 
measurements of qn n k- For that, we need, at each temperature, two copies of the system, so we actually run 2 sets of 
Nt replicas and perform the global moves independently in each of these two sets. 

For the three-dimensional model we stored configurations for sizes L = 4, 5, 6 and 8 at T = 0.20, 0.50 and 2.0, which 
are to be compared with [^0| T c ~ 0.95. We also stored size L = 12 configurations at T = 0.50. The parameters of 
the simulations are shown in Table ffl. The highest temperature was 2.0 and lowest 0.2 except for L = 12 where the 
lowest temperature was 0.5. 

We generated randomly chosen interactions, Jy, with a Gaussian distribution with zero mean and standard deviation 
unity. For each size, temperature and bond configuration (sample) we saved 500 spin configurations. These, together 
with the 500 obtained from them by spin reversal, constitute our ensemble of M = 1000 spin configurations, generated 
for each sample. 

For the four-dimensional model we stored configurations for sizes L = 3, 4 and 5 at T = 0.2, 0.8 and 2.6, compared 
with |34|] T c ~ 1.80. The highest temperature was 2.6 and the lowest 0.2. 500 spin configurations were saved for each 
sample. The other parameters of the simulations are also shown in Table [l| 

We are confident, based on the equilibration test used p3], that the spin configurations we generate are in thermal 
equilibrium. However, it is interesting to ask whether there are significant correlations between them. Our results 
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TABLE I. Parameters of the simulations in D = 3 and 4 dimensions. N sa .mp is the number of samples (i.e. sets of bonds), 
fiequil is the number of sweeps for equilibration and n meas is the number of sweeps for measurements for each of the 2Nt replicas 
for a single sample. Nt is the number of temperatures used in the parallel tempering method. 



do not require that correlations be absent, but the clustering method does require that a substantial number of 
independent configurations are generated for each sample. 

For each set of bonds (and temperature) we store 500 spin configurations, 250 for each replica, so the number of 
sweeps between measurements, t mcas , is given by t mcas = n mcas /250 where n mcas is given in Table | We will denote 
by "time" , t, the number of Monte Carlo sweeps. A quantity which tests for correlations is the time-dependent 
Edwards- Anderson order parameter [g(t)]j = [(Si(t )Si(t + t))]j, where (• • •) indicates a thermal average. This is 
estimated from our spin configurations according to 



W)h = 
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where we have averaged over N to values for the initial time to as well as over spins and bond configurations. Clearly 
[?(0)] j = 1 and j — > for times sufficiently long that there are no correlations. 

In Fig. |^ we show data for [g(£)]j in D = 3 for L = S,T = 0.2. We see that the correlation is very small even 
for t/t meas = 1 (i.e. between the configurations of neighboring measurements). The same is true for smaller sizes 
and higher temperatures. For L — 12, T — 0.5, shown in the inset to Fig. ^, the correlations are larger, about 0.24 
for t/t mcaa — 1, and then decrease, though less fast than exponentially. Thus, for L = 12, correlations will decrease, 
somewhat, the effective number of independent spin configurations. However, we feel that this is not crucial since we 
do not use the D = 3, L = 12 data for the clustering analysis, and only present it in one place, Figure [| 

In D = 4, for L = 3 and 5, the strength of the correlations at T = 0.2 is small, comparable to, or less than that for 
D = 3, L = 8, T = 0.2. For L = 4, the correlation is intermediate between the results shown in D = 3 for L = 8 and 
12. 



III. CLUSTERING METHODOLOGY 



Clustering is an important technique to perform exploratory data analysis. The aim is to partition data according to 
natural classes present in it. By "natural classes" we mean groups of points that are close to one another and relatively 
far from other points, so that it is natural to assign them together, without using any preconceived information on 
the features according to which the set should be classified. 

The standard definition of the clustering problem |35| is as follows. Partition N given data points (or objects) into 
K groups (i.e. clusters) so that two points that belong to the same group are, in some sense, more similar than two 
that belong to different groups. The i = 1, 2, ...N data points are specified either in terms of their coordinates Xi in 
a D-dimensional space (representing the measured values of D attributes or features) or, alternatively, by means of 
an N x N "distance matrix", whose elements dij measure the dissimilarity of data points i and j. The traditional 
tasks of clustering algorithms are to determine K and to assign each data point to a cluster. 

In the context of the present work we can think of our sample of M spin configurations as the objects to be clustered. 
Each object is represented by an N— component vector = (Sf , Sq,—, wnere = ±1 is the value taken by 
spin i in state /i. An alternative view, which we also use, is to consider the N spins as the objects to be clustered. 

Our first aim in this work was to look for a hierarchical structure of the states of a spin glass. Hence we wanted 
to find a hierarchy of partitions, where each partition is a refinement of the previous partition. This purpose calls for 
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FIG. 2. The main part of the figure shows the correlation between spin configurations, [q(t)]j, defined in Eq. (^|) of the text, 
in D — 3 for L — 8, T — 0.2. The horizontal axis represents the number of Monte Carlo sweeps between the two configurations 
in units of the number of sweeps between individual measurements, £ mcas . For comparison, for each set of spins ("replica"), a 
total of 250 configurations are generated. The inset shows results for D — 3, L — 12, T = 0.5, which indicate that correlations 
between spin configurations are significantly larger than for L = 8. 



using a hierarchical clustering algorithm. The output of such an algorithm is a tree of clusters, called a dendrogram. 
Each node in the tree corresponds to a cluster. The splitting of a cluster represents its partition into sub-clusters. 
The trunk is the single "cluster" that contains all the objects, representing the crudest partition; at the other extreme 
each leaf is a cluster of a single object, representing the finest partition. 

There are many clustering algorithms thatproduce such a hierarchical partitioning of any data set. We tried two 
algorithms; a recently introduced one, SPC |3q| , which uses the physics of granular fcrromagnets to identify clusters, 
and a graph-based algorithm proposed by Ward. In the present problem the state clusters are nearly always compact 
(i.e. consist of a high density of points concentrated in a relatively small volume), and the same holds for spin clusters. 
Therefore an algorithm that identifies compact clusters easily is most suitable for our needs and Ward's algorithm 
is designed to find such clusters. Furthermore, SPC is a "short-range" algorithm 37 1, in the sense that it couples 
directly only points within a characteristic length scale. If this scale is tuned by the distances inside valleys, which 
are much smaller than the distance between them, SPC identifies the valleys as different clusters, but may miss their 
relative hierarchical structure. 

Ward's algorithm ]|5[ is agglomerative, works its way up from the leaves to the trunk, by fusing two clusters at each 
step. It begins with an initial partition to i = 1, 2, N clusters, with a single data point in each. One calculates the 
distance Dij between every pair of points one may use, for example, the Euclidean definition of distance, or (for 
binary valued coordinates) the Hamming distance. 

At each step that pair of clusters, a, /3, which are separated by the shortest effective distance p a p from each other, 
are identified and fused to form a new cluster a' = a U f3. The process stops when there is only one cluster, that 
contains all points. 

Initially each data point i = 1, 2, ...N constitutes a cluster and hence the distance pij between two such "clusters" is 
the original distance between points i and j. For subsequent steps, however, one must define an effective distance 
Pa/3, between any two clusters a and (3. This distance is defined by the following update rule: if at a particular 
step we fuse two clusters, a and /?, to form a new cluster a', we calculate the effective distances pL a i, between every 
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FIG. 3. (a) The dendrogram obtained by clustering the M 
at T — 0.2. The vertical axis describes the value of r, defined in Section 
input to Ward's algorithm. Darker shades correspond to smaller distances. The states are ordered according to their position 
on the dendrogram (a), (c), (d) The same as in (a), (b), for the same realization {J}, but for an ensemble of states obtained 
at at T — 0.5. (e), (f) The same as in (a), (b), for the same realization, at T — 2.0, which is greater than T c ~ 0.95. Note 
that this dendrogram is not symmetric; almost all the distances are close to 0.5, so at each stage of the algorithm there were 
several possible partitions that gave minimal value to S. In the implementation we used, the algorithm chose a non-symmetric 
partition. 



unchanged cluster, 7 ^ a, /3, and the new a', according to 



Pa'-y = " Paj H T T PPl j : Pa/3 , W 

n a + np + n~ ( n a + np + n 7 n a + np + n 7 



where n x is the number of data points in cluster x. Distances between unfused clusters remain the same. Note that 
p' a ,~ > p a p and p' lS > p a p for every two clusters 7, S. Hence after every fusion step the minimal distance between 
clusters increases. 

Whenever two clusters are fused, the quantity 



s = J2^ ( 7 ) 



where a a is the sum of squared distances over all pairs of points in cluster a, 

<r a =Yl D » 2 - ( 8 ) 

i,j£a 

increases. It can be shown |$5| that Ward's fusion and distance update rules ensure that at each fusion step this 
increase is minimal. 

We associate a value r with each cluster a', where r(a') = p a p is the effective distance between the two clusters 
that were fused to form a' . For the initial single-point clusters we set r = 0. r(a) is related to a a , the sum of squared 
distances within cluster a. Clusters formed earlier have lower r values, and their a a is smaller. 

The result of the algorithm is a dendrogram, or tree, as in Fig. ||(a). The leaves at the bottom represent the 
individual data points; they are ordered on the horizontal axis in a way that reflects their proximity and hierarchical 
assignment p8[ . The small boxes at the nodes represent clusters. The vertical location of cluster a is its r value, and 
is thus related to its a. When two relatively tight and well-separated clusters are fused, the r value of the resulting 
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cluster is much higher than those of the two constituents. Hence the length of the branch above cluster a provides a 
measure of its relative o~ a ; long branches identify clear, tight clusters. 

Every clustering algorithm is designed to work well for data that satisfy some (usually implicit) assumptions. 
When the actual distribution of the data points deviates from these assumptions, the algorithm may produce some 
"unnatural" partitions. For Ward's algorithm one has to look out for two potential problems. 

The first problem arises from the implicit assumption that minimizing S, the variance within clusters, leads to 
"natural" partitions. This is not the case when, for example, the data consists of a set of points C whose natural 
partition is into two clusters C\ and Ci with very different sizes. We encountered this problem only for the classification 
of very small groups of states, and therefore it has very little statistical effect on our results. 

The second, and seemingly more serious concern is the fact that like every agglomerative algorithm, Ward's algorithm 
will generate a tree-like structure when applied to any set of data. In fact, it is fairly easy to identify when the 
dendrogram and the corresponding partitions do correspond to real hierarchical structure, and when is it an artifact of 
the clustering algorithm used. We used three indicators for the "naturalness" of our state clusters: direct observations 
of (1) the dendrograms and (2) the distance matrices, as well as (3) a quantitative measurement of the sizes of our 
clusters, which are significantly smaller than the distance between clusters. These points are demonstrated in Sec. 



IV; for a detailed discussion see 



IV. STATE SPACE STRUCTURE 

For a particular (randomly chosen) set of bonds {J} of the system we generate, as discussed in Sec. H a sample of 
500 states, which constitute an equilibrium ensemble at a temperature T. Next, we add to this ensemble the set of 500 
states obtained from the original set by spin reversal. Clearly the new ensemble of M = 1000 states also corresponds 
to thermal equilibrium p9|] at T. We now address the following question: 

Do the M states of the equilibrium ensemble cover the 2 N points of state space or a part of it uniformly, 
or is there some underlying hierarchical organization? 

As it turns out, the answer depends on T; whereas above T c the M states do not exhibit any apparent structure, 
below T c a very pronounced hierarchical organization is seen. To uncover this organization we use the clustering 
methodology of the previous Section, treating the M states of our ensemble as the data points to be clustered. 

We describe here analysis of a single realization of the randomness, in order to help the reader perceive the qualitative 
nature of the results (see Figs. [| and ||), and to define the observables that we measure. These observables were 
measured for each of the different realizations, and the distributions of their values were determined; the average and 
width of these distributions are also presented. This data demonstrate that the results described in this section for a 
single sample are typical and seen in many samples. 

In order to cluster the states, each state fi is represented as an iV-component vector S M = (S^, S^r), where 
= ±1 is the value taken by spin i in state /i. The complete data set can be represented as an N x M data 
matrix, whose columns are the vectors S M . For the set of M = 1000 states, obtained at T = 0.2 for a particular bond 
realization of an N — 8 3 spin system, the data matrix is presented in Fig ||(a). Pixel of this figure represents 

the sign of spin i in state y,; a black entry corresponds to +1 and white to — 1. The spins appear in lexicographic 
order and the states in the random order generated by the simulation. As can be seen, the matrix appears fairly 
random, with no easily discernible structure; nevertheless, there is a clear organization of these M states into tight 
clusters. For the particular realization and ensemble of states presented here, these clusters of states can be seen by 
direct observation of the M = 1000 data-points S M , once one overcomes the hurdle of directly viewing a cloud of 1000 
points in a N = 512 dimensional space. 

A trivial way of visualizing points that lie in a high dimensional space is to project them onto a low (i.e. two 
or three) dimensional subspace. In order to reveal the underlying structure, it is important to choose with care the 
subspace onto which one projects. A widely used method to choose this subspace is that of principal component 
analysis (PC A) fib]] . One constructs the N x N covariance matrix of the M points, 

M 

^■^E^" ( 9 ) 

where 

5S? = {S? - mOto (10) 
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with rrii the average of the M variables Slf and of their variance. For our case rrii = and <jj = 1 for all i, and hence 
the covariance matrix is the spin correlation matrix, i. e. 



1 M 



(ii) 



The eigenvectors e$ of this matrix are the principal directions or components of the variation in the data. They are 
ordered according to the size of the corresponding eigenvalues, with the largest coming first. 

In Fig. H we present the projections of our T — 0.2 ensemble of M = 1000 states on the first two and three principal 
components. Even though projection of N = 8 3 dimensional data onto three and two dimensions involves a major 
loss of information, the cluster structure of the states is still clearly evident. In Fig. || (a) projection onto the largest 
eigenvector, ei, is represented by the horizontal axis, and on the second largest, e2, by the vertical. It is interesting 
to note that the two largest state clusters, C\ and C±, project mostly onto ei and the second largest pair, Ci and C2 
onto e 2 . Fig. |5| (b) indicates that the next sized variation, due to splitting of C2 into two subgroups, is captured by 
63. The scale of the projections can be understood by the following argument: if the (normalized) eigenvector ei is 
parallel to a typical vector from Ci, then, since normalization of ei involves a factor of 1/y/N, the maximum possible 
projection is y/N ss 22.6. Hence the projections shown in Fig. ^| are quite large, i.e. close to the maximum possible 
value. 

Next we obtain a systematic quantitative measure of the hierarchical structure of state space by performing a cluster 
analysis of the M points. The choice of the particular clustering algorithm used was dictated by our idea of the state 
space structure, obtained from PCA and from the picture described in the Introduction and summarized in Fig. [l]. 
Our aim is to find a hierarchy of partitions into compact clusters. That is, we would like states that belong to the 
same cluster to be closer to each other than to states in different clusters. Ward's algorithm, described in Section p| , 
is tailored to perform this task for the kind of data distribution that we have in state space. 

To start, we defined the M x M distance matrix D between the states /i, v by 

B. v = , (12) 

where q^ LV is the state overlap defined by Eq. (^5j). Next, we clustered the states using the distance matrix D^ u as input 
to Ward's algorithm (see Eq. @). The algorithm results in a dendrogram, as shown in Figs. ||(a,c,e), for a sample at 
T = 0.2, 0.5 and 2.0, in three dimensions. The leaves, which represent the states, are ordered on the horizontal axis 
according to the order imposed by the dendrogram jj8). The nodes represents the clusters. The vertical location of 
each cluster corresponds to its r value, and is thus related to the variance within it. 

For T = 0.2 and 0.5, which are below T c w 0.95 ^0|, we found clear partitions in the two highest levels of the 
dendrogram, as presented in Figs. ||a,c). At the highest level the states are partitioned into C and C. At the next 
level, C is broken into two sub-clusters, which we denote as C\ and C 2 . For this specific sample the cluster C 2 breaks 
further into two sub-clusters, which are clearly seen in Fig. |H| as well. 

To gain insight into the manner in which similar states are grouped together, and to actually "look into the spin- 
glass" at the microscopic level, we present in Fig. [|(b) the same data matrix as shown in Fig. f|(a), but with the 
states again reordered according to the dendrogram of Fig. I (a). That is, to get Fig. | (b), the columns of Fig. |(a) 
have been permuted according to their position in the dendrogram. The clear central vertical dividing line separates 
C from C. In addition to the central dividing line, another vertical line is also clearly visible - it separates the states 
that belong to the larger cluster C\ from the smaller one, C 2 . 

We now demonstrate that the state clusters we found are indeed "correct" and "natural". First, we checked that 
the situation of merging two clusters of very different sizes occurs very rarely. 

We showed that our partitions are "natural" and not an artefact of the algorithm (which produces a tree for any 
data), in three ways: 

1. Note that direct observation of the dendrograms clearly distinguishes between the different situations above and 
below T c . At T = 0.2,0.5 (< T c ) the relative r values of the state clusters C,Ci an d Ci ~ measured by the 
length of the branch above each cluster - are high. A long branch indicates that the size of the cluster is much 
smaller than the distance between it and its "brother" , which indicates that the partition into these two groups 
is natural. In comparison, in the dendrogram obtained at T — 2.0 (> T c ), the relative r values are much smaller 
than at T = 0.2,0.5. 

2. The genuinely hierarchical structure at T = 0.2,0.5 is also evident from the states' distance matrix, as shown 
in Figs. ||(b,d). This distance matrix was obtained by reordering the states according to the results of the 
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FIG. 4. (a) The original data matrix of 500 x 2 states S M , S£ = ±1, with black/white representing +/— . This 3D sample 
was generated for a realization of size 8 3 at T = 0.2 (the same one as in Fig. ^|). The spins are in lexicographic order, (b) The 
same matrix, with the states ordered according to the dendrogram in Fig. til (c) The matrix in (b) with, in addition, the spins 
ordered according to the spin dendrogram T> in Fig. O. 



cluster analysis, i.e. according to the order of the leaves of the corresponding dendrogram. When the states 
are randomly ordered (like in Fig. ^(a)), the resulting distance matrix is a homogeneous greyish square, like 
that of Fig. ||(f). The difference between this and Figs. |](b,d) is striking: the distance matrices within clusters 
C\ and C2 appear as dark squares (representing shorter distances) along the diagonal. The distances between 
clusters are represented by fairly uniform, lighter colored rectangles. In comparison, for T — 2.0 there is no real 
hierarchical organization of the states, and reordering them according to the dendrogram does not generate any 
ordered appearance of the distance matrix. 

3. We measured the average distance between pairs of states that belong to each of the clusters C, C\ and C2. The 
average D(C) and the width w(C) of the distribution of distances within C are 

£>(<?) = TcfpE^ee A- ! (13) 
HC) = E»,„ec D ^ 2 - D (C?) V2 , (14) 
where \x and v refer to individual configurations. The average D(C a ) and the width w(C a ) for a = 1,2 are 
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FIG. 5. Principal Component Analysis of a sample of M = 500 x 2 states of a specific realization of {J} in 3D with N = 8 3 
spins at T — 0.2. Each point represents a state S M . The coordinates are projections on to eigenvectors e*, corresponding to 
the largest eigenvalues of the correlation matrix in Eq. ([n]). We show in (a) projections onto two eigenvectors, corresponding 
to the largest and next-largest eigenvalues of the correlation matrix, shown, respectively, on the horizontal and vertical axes. 
In (b) the three largest eigenvectors are used. The first and second level partitions of the hierarchy are clearly visible and, to 
some extent, the third level also. 



defined in a similar way. The distribution of distances within clusters is to be compared with the distribution 
of distances between points that belong to different clusters. The average D(Ci,C2) and width w(C\,C2) of the 
inter-cluster distance distribution are defined as 



D ( C ^^) = Jc^cT\^ec 1 E,ec 2 D^; (15) 
^Ci,C 2 )= (^£ pGCi ^ eC2 IV 2 -#(Ci,C 2 ) 2 ) 1/2 . (16) 



The clusters C,C are special in that each state \i € C has an inverted state /j, G C, so that S' 1 = — S**. Therefore 
D(C,C) = 1 - D(C) and w{C,C) = w{C). 

A subset of the results is presented in Table O; for all temperatures, system sizes and both dimensions see [p7| . 
We present for each variable x its mean [x]j (averaged over the disorder {J}) and its standard deviation 
Ax = ([x 2 ]j — [xjj 2 ) 1 / 2 . For T = 0.2 and 0.5, which are below T c , the average distances within the clusters 
are of the order of 0.1. D(C,C) is around 0.9, which shows that there is a clear separation between these two 
clusters. D(C\ : C2) is much lower, but is still about two or three times larger than either D(C\) or D(C2)- Note 
that the width of the distance distribution within a cluster is of the same order of the mean distance, so in 
general distances will not be much larger than twice the mean distance. At T = 2.0 (> T c ) the distances 
within and between clusters are almost equal and the differences are only due to statistical fluctuations, again 
indicating absence of natural structure, as we claimed on the basis of direct observation. 

Measurement of some of the quantities listed above allows us to investigate the extent to which the state space 
structure of short-range spin glasses, as reflected by the data in Table ||, is compatible with RSB. In the RSB [|-||] 
framework, the overlap between any pair of valleys (which correspond to pure states in the usual interpretation of 
RSB) from two different clusters that appear at the same level of the hierarchy is constant. It seems natural to 
associate the pure state clusters of RSB to our state clusters, e.g. C\ and Ci- In this association, each state cluster 
contains states that belong to different "pure states" . If the overlap between pure states of the two clusters is constant 
as in RSB, this should hold also for the overlap between each pair of states \i G C% an d v € C2, since the width of the 
overlap distribution inside a pure state approaches zero. In this case, all entries of the sub-matrix for /x G C\ and 
v G C2 would be equal, so the width 1012 = [wiCijC^)] j should vanish as L — > 00. To test whether this is the case, 
we present in Fig. ^ the values of w\2 = [w(Ci, C2)],/ vs the system size L for T = 0.2 and D = 3. The error bars 
represent the statistical error (obtained by dividing the standard deviations, given in Table ||, by y/N sa , mp — 1). We 
tried fits of the form 

wn = Woo + BL-y , (17) 

with B and y as fit parameters. The overall best fit was for Woo = 0.0205, B = 0.58 and y — 3.36, which gives a very 
small x 2 °f 0.036. This is shown by the solid line in Fig. H. We also tried the best fit assuming that w^, = 0, which 
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FIG. 6. A log-log plot of Wi2 against L for T = 0.2 and D = 3. The solid line is the best least squares fit to Eq. (|l7|), while 
the dashed line is the best fit with the additional assumption that Woo — 0. . 

has fit parameters B — 0.039 and y — 0.30, and is shown by the dashed line in the figure. This has a x 2 of 1.41 which 
is much larger than the best fit with Woo ^ 0, but still acceptable. Hence even though our data suggests that iWoc 7^ 0, 
the possibility that = 0, which corresponds to RSB, cannot be ruled out. 

V. CORRELATED DOMAINS IN SPIN SPACE 
A. Identifying the spin domains 

Accordin g to our picture, splitting of a cluster at level a in the states hierarchy is induced by a macroscopic 
contiguous EjJ spin domain Q a . The size and shape of this domain determines the energy barrier separating two 
state clusters that were "born" at this level. In this subsection we describe how we identify from our data the two 
correlated domains Q\ and O2, which determine the two highest levels of the states hierarchy, and also discuss whether 
they remain macroscopic at large L. Domains that emerge at the next level, Q3 and C/3, are also discussed briefly 

Since the spins in such domains flip "collectively", they are highly correlated. The standard definition of the 
correlation Cy of spins i and j is 

en = (S t S 3 ) = I Yl S i S i exp[-/3W(S)] , (18) 
s 

where (...) stands for the thermodynamic average for a particular realization of the disorder, and 2 is the partition 
function at T. Using our equilibrium ensemble of states {S M }, we evaluate 

^4w- (19) 

The correlation in itself is unimportant for spin glasses since it is gauge dependent and its average [cy] j over all the 
realizations of the disorder { J} vanishes. The relevant measure of correlations in a spin glass is the square, cy . If 
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two spins are independent of each other over the equilibrium ensemble of states, we have Cy 2 = 0. On the other hand, 
for a pair of fully correlated spins we have Cij = 1; the two spins are either aligned or anti-aligned in all states. 
To proceed, it is convenient to define, quite generally, Q^ v as the set of spins whose sign is different in [i and v, i.e. 

G, v = {i\ + S» } . (20) 

We expect the largest domain, Gi, to be in one orientation in the states of C and in the reversed one in the states of 
C. To identify the spins that indeed behave this way, we took all (M/2) 2 pairs of states /i 6 C and v 6 C and, for 
each pair, determined G^u- Ideally all the spins of Gi always flip together and maintain their relative orientation; if 
so, the set of spins G^v for all pairs of states fi and v would always include Q\. However, at finite T we must allow 
for excitations of the order of J. So, even if a spin is highly correlated with the other spins of Gi, it might lose its 
relative orientation in a few of the M states of the sample. In order not to "miss" such spins, we use a soft criterion 
when we determine whether a spin is a member of Q\. We define a threshold 9 and define Gi{9) as the set of spins i 
which are members of G^, i.e. for which S^S\ — —1, for at least a fraction 9 of the pairs of states fi € C and v € C. 
This can be written as 



since the terms in the normalized sum where S^S\ = 1 must, by definition, sum up to less than 1 — 9 and the sum 
of the terms with S^S" = — 1 must be less than —6. We define our spin domain Gi(9) as the largest contiguous part 
of Gi{9). For large enough 9 we found that for most realizations {J}, below T c the sites of Gi(9) are contiguous and 
hence it is identical to Gi(9) (for detailed values of the ratio |<? a |/|£7i|, its mean over realizations and its standard 
deviation, see fl37f ). The next spin domain G2(9) is defined in the same manner, on the basis of pairs of states fj, € C\ 
and € C2. 

The above definition sets a lower bound on the correlation of spins within the domain. Consider two spins i,j S 
Gi(9). By definition, 

^ 2 -^E^^- ( 22 ) 

Now the number of states in C and C are both equal to M/2. In addition, for a given v, we can replace fi by its 
inverse p, and the product of the four spins doesn't change. Hence we get the same contribution from fi £E C as \x G C. 
As result we have 

Cy2 = |C||C| ^ ^ S '' S '' S '' S j ■ ( 23 ) 

Now S^Si will be —1 for a fraction of the states fj, and v which is greater than 9 and +1 for a fraction less than 1—6, and 
similarly for SjSj. Hence S^S" and S^S'j will have the same sign with probability greater than 1 — 2(1 — 0) = 29—1. 
Consequently, for i,j£ Gi(9), we have 

Cij 2 > 29 - 1 - [1 - (26» - 1)] =46»-3. (24) 

The same constraint holds also for G2, with the sums taken over the states in clusters C% and Ci- 

Since we introduced an arbitrary parameter 9 into the definition of our spin clusters, it is important to consider the 
extent to which the value of 9 affects their identification. As seen in Fig. ^, the sizes of the domains and their average 
correlation, defined below in (B5|), do not change much for 0.6 < 9 < 0.95. For both a = 1,2 we define (arbitrarily) 
Ga = G a (0.95). We do not choose 9 = 1 since, as discussed above, we do not want our results to be affected by small 
thermal fluctuations. In Fig. || we plot the spatial structure of Gi and Gi for a specific realization. For T > T c 
the correlations between each pair of spins are much smaller, and hence this analysis is meaningless. The procedure 
described above results in Gi{9) = G2{9) = for any 9 > 0.5. 

According to our picture these correlated spin domains govern the hierarchical structure of state space. It is 
important to clarify whether these domains survive as the system size L increases. There are two mechanisms by 
which increasing the system size can invalidate our picture: either the domains do not remain macroscopic when L 
increases, or they do remain macroscopic but merge as L —> 00, i.e. the fraction of states in which Gi flips tends 
to zero. We now discuss each of these possibilities in turn. In addition, a simple figurative description of these two 



mechanisms is given Sec. VB 
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1. The domains do not remain macroscopic when L increases. To study the finite size effects of our analysis we 
normalized the domain sizes by the number of spins and plotted the size distributions of the two domains for 
different system sizes, in D = 3 (see Fig. ^ and in D = 4 (Fig. 10), at two temperatures in both dimensions. 
The number of bond realizations, iV samp , from which these distributions were obtained for various system sizes 
i, at both D = 3,4, are given in Table B. For T = 0.2 in both dimensions, and at T = 0.5 for D = 3 the 
distributions seem to converge even for the small system sizes we use. We conclude with high certainty that at 
T = 0.2 for D — 3, 4 and at T = 0.5 for D — 3 the domain sizes \ Q a \ are proportional to L D for both a = 1, 2. 



The mean and width of these distribution are presented in Table III . The width of the distributions does not 
vanish, so the sizes of the domains are non self-averaging quantities. On the other hand, for T = 0.8 in D = 4 
we cannot determine conclusively whether the domain sizes do or do not remain proportional to N = L D as L 
increases. 

2. Q\ and Gi may remain macroscopic but merge as L — > oo. If this occurs, we end up with a single domain and 
there will be no hierarchical structure in state space. To check that this does not happen we calculated the 
average correlation C12 between spins in Gi and C? 2 j 



1 

Cl2 



01 1 1& 



E E c - 2 • (25) 



If C12 approaches the value 1 as L — > oo, the two domains indeed merge in the thermodynamic limit. In Table HI 
we present, for systems of different sizes and dimensions, the average values of c.\i (averaged over the disorder 
{J}) and the corresponding standard deviations. For T = 0.2, D = 3,4 and for T = 0.5, D = 3 the average 
correlation decreases slightly as the system size increases, although, in D = 3 it seems to converge already for 
L = 8 to a fixed value of ~ 0.5. This means that the spins of Gi and Q2 will not become fully correlated and 
the two domains will stay separate as L increases. 

Interestingly, in D = 4, the correlation for L = 4, 5 is higher at T = 0.8 than at T = 0.2. The reason for this 
is probably that as T increases, small pieces of Gi "fall of" . Since G2 at T — 0.2 is small, one of these pieces, 
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T 


L 


[D(C)]j 


AD(C) 


HC)]j 


Aw(C) 


[D(Ci)]j 


AD(d) 


[w{d)]j 


Atu(&) 


0.2 


4 


0.045 


0.049 


0.055 


0.052 


0.015 


0.017 


0.019 


0.018 




5 


0.050 


0.054 


0.056 


0.054 


0.018 


0.018 


0.019 


0.019 




6 


0.053 


0.056 


0.054 


0.053 


0.021 


0.020 


0.019 


0.019 




8 


0.055 


0.054 


0.052 


0.051 


0.025 


0.020 


0.020 


0.020 


0.5 


8 


0.139 


0.065 


0.084 


0.046 


0.093 


0.038 


0.045 


0.026 




12 


0.151 


0.065 


0.078 


0.046 


0.106 


0.036 


0.041 


0.024 


2.0 8 


0.487 


0.006 


0.053 


0.002 


0.477 


0.009 


0.055 


0.002 




[D(C 2 )]j 


AD(C 2 ) 


MC 2 )]j 


Aw(C 2 ) 


[D(C u C 2 )]j 


AD(Ci,Ca) 


K&,Ca)]j 


Aw(C u C 2 ) 


0.2 


4 


0.025 


0.036 


0.027 


0.034 


0.160 


0.135 


0.026 


0.024 




5 


0.025 


0.032 


0.025 


0.031 


0.169 


0.147 


0.023 


0.020 




6 


0.028 


0.033 


0.026 


0.033 


0.161 


0.141 


0.022 


0.021 




8 


0.030 


0.027 


0.024 


0.026 


0.161 


0.139 


0.021 


0.018 


0.5 


8 


0.112 


0.057 


0.057 


0.037 


0.253 


0.126 


0.053 


0.027 




12 


0.121 


0.048 


0.054 


0.033 


0.263 


0.125 


0.044 


0.023 


2.0 


8 


0.472 


0.009 


0.057 


0.002 


0.499 


0.005 


0.048 


0.003 



TABLE II. The average distances within and between state clusters, and the relations between them, for a subset of the 
D = 3 dimensional systems. For each variable x we present the average over all realizations, [x]j, followed by its standard 
deviation, i.e. Aa; = ([a; 2 ] j — [as] j 2 ) 1 ^ 2 . The statistical error of each mean [x]j is Ax/ */ Af samp ; the number of samples for each 
L, D is given in Table 



D 


T 


L 


[\Si\]j/N 


AIG^/N 


[\S-2\}j/N 


A\g 2 \/N 


[ci 2 ]j 


Acia 


P{Q 2 + 0) 


3 


0.2 


4 


0.70(1) 


0.21 


0.099(4) 


0.087 


0.56(1) 


0.33 


0.856(6) 






5 


0.66(1) 


0.21 


0.105(5) 


0.104 


0.55(1) 


0.33 


0.832(6) 






6 


0.66(1) 


0.20 


0.090(4) 


0.090 


0.52(2) 


0.34 


0.836(6) 






8 


0.64(1) 


0.20 


0.084(5) 


0.094 


0.53(2) 


0.34 


0.833(8) 




0.5 


4 


0.31(1) 


0.21 


0.062(3) 


0.056 


0.49(1) 


0.32 


0.56(1) 






5 


0.26(1) 


0.18 


0.052(2) 


0.043 


0.49(1) 


0.33 


0.57(1) 






6 


0.25(1) 


0.16 


0.046(2) 


0.046 


0.47(1) 


0.33 


0.52(1) 






8 


0.22(1) 


0.15 


0.035(2) 


0.034 


0.47(2) 


0.31 


0.55(1) 






12 


0.24(1) 


0.15 


0.033(2) 


0.035 


0.54(2) 


0.31 


0.56(2) 


4 


0.2 


3 


0.74(1) 


0.19 


0.107(5) 


0.105 


0.62(2) 


0.34 


0.840(6) 






4 


0.73(1) 


0.19 


0.083(4) 


0.092 


0.53(2) 


0.34 


0.830(6) 






5 


0.73(1) 


0.19 


0.082(7) 


0.098 


0.51(2) 


0.34 


0.77(1) 




0.8 


3 


0.154(7) 


0.15 


0.036(1) 


0.031 


0.47(1) 


0.31 


0.298(9) 






4 


0.142(6) 


0.12 


0.025(1) 


0.029 


0.54(1) 


0.31 


0.37(1) 






5 


0.139(8) 


0.11 


0.020(2) 


0.025 


0.57(2) 


0.29 


0.38(2) 



TABLE III. The normalized sizes of the domains Q\ and Q 2 , and the average correlation between spins that belong to the 
two domains. The last two parameters are taken for realizations {J} where Q 2 does not vanish. The probability for Q 2 not to 
vanish is also presented. For each quantity x the table contains [x]j, its average over JV sam p realizations of the disorder {J} 
and the width of the distribution Aa; = \/[a; 2 ]j — \x\j ■ Next to each [x]j we show its statistical error (in parentheses). 
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FIG. 8. The spin domains Q\ and Q2, as found in the realization of Fig. ^ Note that we use periodic boundary conditions, 
so the domains are connected through the boundaries. No spin is shared by Q\ and C/2- 




which is larger than G2, assume the role of G2 at T = 0.8. Since this piece was part of Gi at T = 0.2, we expect 
its correlation, with what remains of Gi at T = 0.8, to be relatively high. Extrapolating from L = 3, 4, 5 is not 
useful, but we still believe that the correlation does not approach 1 as L — ► 00. 

We also attempted to identify G3 and G'3, the spin domains associated with the third level of the state hierarchy 
(see below). G3 is the cluster which is associated with splitting C\ into its two descendents on the dendrogram, C\ a 
and Cib. The domain G'3 plays the same role in C 2 - Since by our notation \C\\ > \C<z\ we expected that in order to 
have a larger number of states, the spin correlations will be lower when measured over C± than over C%. As a result we 
expect 1 03 1 < |0g|. Due to the small sizes of the systems we study, we cannot be sure if the sets of spins we identify 
as G3 and G' 3 indeed play the role we attribute to them, or are just a microscopic noise and, therefore, only a finite 



size effect. The results are given in Table [TV . We see that the normalized sizes of both domains decrease with the 
system size, perhaps due to finite size effects. We also measure the average correlation c(G3, Gi U G2), of G3 with the 
largest domain correlated over C\, which includes Gi U G2 (this domain has a fixed orientation over the states of C\). 
This correlation is defined as 

^ giUQ2) = w^km E (26) 



In Table IV we see that the values of c(G3,Gi U G2) decrease as L increases; hence if G3 survives as a macroscopic 



1G 




cluster at large L, we expect it to remain distinct from the union of the two larger domains. 



B. Spin space structure 



So far we have obtained the spin domains using the results of the state space analysis. However, the existence of 
these domains can also be observed directly in spin space, i.e. without utilizing information about the previously 
identified hierarchical structure of state space, as we now demonstrate. 



As described in Sec. IV, the equilibrium ensemble of states, obtained for each realization, is represented by an 



N x M data matrix {Sf } (e.g. Fig [|(a)). In Sec. IV we treated each of the M states, represented by a column of this 
matrix, as a "data point" whose coordinates are the components of this TV-dimensional vector. Now we view each of 
the N spins of the system as a data point, represented by a row of the same matrix. Each of these data points is a 
vector in an M-dimensional space. 

The distance on the set of spins should be defined according to the nature of the clusters we are interested in. At 
this case, we expect highly correlated spins to be in the same cluster, and spins with low correlation to be in different 
clusters. Thus, we define the distance between a pair of data points i and j as 

dij = 1 - dj 2 . (27) 

This N x N distance matrix serves as the input for clustering the spins, using Ward's algorithm. The dendrogram 
V, obtained when the data of Fig. |](a) are clustered, is presented in Fig. |ll|(a). The correlated spin clusters are 
represented by boxes in the dendrogram - let us denote them by g a . When the spins are reordered according to the 
dendrogram, their distance matrix, shown in |ll|(b), clearly exhibits a non-trivial structure. There are large, highly 
correlated spin clusters on the lower levels of the dendrogram. 

In order to "see" the manner in which the spins are ordered, we return to the data matrix of Fig. ^(a). We obtained 
Fig. |](b) from (a) by reordering the columns according to the state dendrogram in Fig. ^. If we now reorder the 
rows of Fig. |^(b) according to the spin dendrogram V in Fig. [ll], we get Fig. [|(c), which is redrawn as Fig. [l^ with 
labeling of the largest state clusters and spin domains. The cluster structure of the spins can be clearly be seen in 
Fig. [[J]. Spins in Qi clearly have the same orientation in the states of C but are inverted in the corresponding states 
of C. Spins in Q2 have opposite orientations in C\ and C2 and are inverted in the corresponding states of C\ and Ci- 
One can also see that spins in domain Q' z separate C2 into two sub-clusters. As to G3, we point in Fig. ^ to a few (3 
- 4) spins, which have the same sign in all states of C2 but change sign in C±. 
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FIG. 11. (a) The spin dendrogram T> for the data of Fig. ^(a) produced by Ward's algorithm, (b) The spin distance matrix 
d of this realization realization. The spins are ordered according to their clusters in T>. Darker shades correspond to smaller 
distances and higher correlations, (c), (d) The same as in (a), (b), for the same realization at T = 0.5. (e), (f) The same as 
in (a), (b), for the same realization at T = 2.0. The j/-axis is rescaled to show the dendrogram, which clearly differ from the 
dendrograms in (a) and (c). 



These data were obtained at T = 0.2 (< T c ). Above T c the correlation between any two spins is low, and there is no 
cluster structure, as evident from Fig. |ll|(e,f). The relative r values of this dendrogram are much smaller then those 
of the dendrogra ms in Figs. |ll](a,c), and the reordered distance matrix is structure-less. If the domains Q a (that were 



identified in Sec. V A on the basis of the state hierarchy) are not an artifact of our analysis, they should be clearly 
identifiable in spin space, and appear as clusters in the spin dendrogram T>. To check this, for each realization we 
compared every spin cluster g a , that appears in the corresponding spin dendrogram D, to every spin domain Q a that 
was previously found for that realization. The spin cluster g a that was found to be most similar to Q a was identified 
and denoted by g a . We used the similarity measure 

o{g a ,ya) = -, — ,, r (,2a) 

\9a\ + \Ga\ 

which represents the fraction of shared spins by the "physical spin domain" Q a and the spin cluster g a . For most 
realizations we have (at low T) g a = Q a for both a = 1, 2; and when these groups are not precisely equal, they differ 
by only a few spins (see |5?J for full details). 

Fig. 03 also provides a convenient, simple "geomet rical" interpretation of the two tests for the survival of our 



picture in the large L limit that we discussed in Sec. VA. Observe the rectangular region corresponding to spin 



domain Q2 and state cluster C2. Validity of our picture relies on "survival" of this rectangle as we take the L — > 00 
limit. The first test we performed checked whether its vertical side, \Q 2 \ stays finite. If this condition is not satisfied, 
the relative area of our rectangle goes to zero; a non-vanishing limiting \Q%\ does not, however, guarantee that the 
rectangle stays finite; it may disappear if its horizontal dimension shrinks to zero when L — > 00. The second test, 
showing that the correlation cyi does not approach 1, ensures that this does not happen either. 

Overall, Fig. n% summarizes in a convenient pictorial way our picture of the spin glass state in short range systems. 



C. Spin domains and states hierarchy 



Now that the spin domains have been well defined, we can examine the manner in which they govern the hierarchical 
partitioning of state space. Each state cluster at level a of the hierarchy can now be identified with one of two possible 
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FIG. 12. A redrawing of the ordered data matrix of Fig. tac), in order to highlight the state clusters and spin domains 
discussed in the text. It is for a 3D realization of size N = 8 at T — 0.2. The columns represent the states /x and the rows 
represent the spins i, = ±1, with black/white representing +/—. The states are ordered according to the dendrogram in 
Fig. ^, and the spins ordered according to the spin dendrogram D in Fig. The state clusters and the spin domains are 
marked (see text). 



configurations of the particular spin domain Q a . We denote these two configurations as -f|- a and JJ-a- Note that we have 
avoided the notation +/— for the states of the spin domains, since in each state some of the spins have the + sign 
and others — . For example, in the first level partition Q\ has a certain characteristic configuration, ffi, over all the 
states in C, whereas over all the states of C it is in the spin inverted configuration 4J- i - The value [lY i ] i , taken by spin 
* G Gi in the configuration ffi, is defined by 



[1h]<=sign VSf ■ ( 29 ) 




Our definition of Qi, using Eq. ( |2l"| ) with 9 — 0.95, guarantees that the argument of the sign function in the above 
expression does not vanish. Hence, stating that Q\ takes configuration \\ in a certain state \x implies that 



ieSi 



S?Mi > . (30) 



The configuration assumed by Q\ in any state /i determines that \i is assigned to C if Q\ is in configuration fti, or to 
C if Gi is in configuration 4J-i . 

The spin domain Qi determines, in a similar way, the partition of C into C\ and C2 (and the partition of C into 
C\ and C2). G2 is in configuration ^2 in states C\ and C2, and in 4J-2 in states C2 and C\ (see Fig. [l] for a schematic 
illustration of this point). 

Each spin domain Q a defines a partition of the states, at level a, into two groups - one in which Q a is in the ff" a 
configuration and the other with -|| a . Picking a pair of states fj, and v, one from each group, the set of spins Qn V , that 
are flipped in the transition between them, will always include Q a [Q. Thus, the distance D^ v = \Q^\/N between 
two such states will almost always be larger then |C* a |/iV. 

By our definition of Q a , the probability that a large part of its spins will lose their relative orientation is small. 
Considering local dynamics, the time it will take Q a to flip is exponential in its size. If Q a is macroscopic (as we have 
shown for a = 1,2) it may be associated with a macroscopic free energy barrier. In an infinite system it will take an 
infinite time to flip, thus inducing a separation of the phase space into two ergodic sub-spaces (or valleys). 
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The clear hierarchical organization of the state clusters suggest that the average distance ( |l5| ) between state clusters 
formed at a high level of the hierarchy is significantly larger then the average distance between clusters formed at a 
lower level. Indeed, we show in Table || that in general D(C, C) ^> D(C\ , C2). We relate this characteristic of the state 
structure to the large variability of the spin domain sizes \Ga\- Indeed, we have seen that typically \Gi\ > 8|(?2| for 
T = 0.2, D = 3,4. 

Now we have a complete picture, supported by our numerical findings, of a hierarchy of state clusters. The valleys 
are the leaves of this hierarchy j|3) . At each level a of this hierarchy the partition of the states is refined according 
to the orientation of macroscopic spin domains Q a . At different nodes of a certain level of the hierarchy there might 
be different correlated domains that determine their partition. Take, for example, the states in C\ (where Gi is in 
configuration \\ and G2 is in configuration ^2). Over these states the largest unlocked jf4| correlated domain is 
G3 = £^3 ^ IT2 ) ■ The two possible configurations of Q 3 inside C\ may be denoted as tf 3 (1T1 , IY2) and JJ.3 (iTi, IV2) • Over 
the states of C2 we expect to find a different unlocked correlated domain GJt = G^Oli, -4J-2 ) - We calculated the part of 
each domain which is included in the other. The results are given in Table IV. We see that G3 and G3 share in general 
less than a fifth of their spins. 

Note that in the ideal case (corresponding to 8 = 1), a spin domain Ga{1li, f|"2, •■■ ilk), that appears at a particular 
level of the hierarchy, cannot share spins with the higher level domains 6=1,2, k, whose orientation is fixed while 
G a flips. For 6 = 0.95 such sharing was also practically excluded. On the other hand, two domains such as G3 and G'3 
can have shared spins, namely those that are free to flip in both the (ffi, 1Y2 ) and (f|~i, -IJ-2) situations. 

Going all the way down the states hierarchy, we find that each valley can be characterized by a specific list of 
domain configurations, e.g. {fti, J| 2 , JJ-3 Oh, h), fN (1hi h, U), ■■■}■ 
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FIG. 13. The two principal components of the 512 spins of three realizations A, B and C (see text) in 3-D. Each point 
represents a spin i and its coordinates are the projections of S; — (Sj , Sf , ■ ■ • , Sf' 1 ) on to the two largest eigenvectors of the 
matrix R in Eq. ([n|). The analysis is carried over (a) all states; (b) the states of C\\ and (c) the states of C2. The spins of Q\ 
are marked by O; of Qi by A; of Qz by Oi an d °f ^3 by x . Spins that belong to both <?3 and G' 3 axe marked by (^). Spins that 
do not belong to any of these domains are marked with dots. The lower half of the plane is projected onto the upper half using 
(x, y) — * (—x, — y). Spins in a correlated domain usually have the same values for the two principal components, and they fall 
on top of each other on the plot. Therefore, in most plots, a correlated domain seems to be represented by a single marker. 



An additional insight is obtained from a PCA of the spins, which is to be distinguished from the PCA of the states 
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in Fig. To perform the PC A of the spins we form the covariance matrix 

1 N 

»=i 

which is analogous to Eq. (|ll]), and project the two largest eigenvectors of R on to the spin configurations Si for each 
site i. 

The results for three realizations, labeled A, B and C, are shown in Fig. [HJ Each data point represents one spin. 
Realization A is the one whose data matrix is shown in in Figs. |(c) and Fig. || 

In the upper left frame of Fig. [l3] we see the results of the PCA analysis of the spins for realization A. We want 
highly correlated spins to be close on the plot. Since a spin Si is fully correlated with its inverse —Si each point (cc, y) 
with y < is projected on the plot to {—x, —y). The spins of Gi are highly correlated with each other and all have 
the same values for the first two principal components of the spin space. Therefore they fall on top of each other, and 
we see only one O marker which represents all of them. The same is true for the spins of Gi, marked by A. As seen 
from Fig. |l2| the spins of Gi are not correlated with the spins of G2 over the M states, and indeed the two domains 
are far from each other on the plot. 

In column (b) of Fig. |l3| we used only the states of C± in the analysis. We can see in Fig. |l2| that over C\ the spins 
of Gi and Gi are correlated, together with some of the spins of G'3, marked by x. In the plot (the middle frame on 
the upper row of Fig. [l3]) we can see that indeed these spins are all plotted at the same coordinates. The spins of £3, 
marked as 0> are highly correlated, but are not correlated with Gi and Gi- Note that the spins of G'3 are separated 
into two different sets, and are not correlated over C± . 

When we perform the analysis using only the states of C2 we get the results presented in column (c) of Fig. In 
the matrix of Fig. [l^ we see that the spins of Gi, Gi and G3 are correlated together over Ci, and indeed they all fall 
on top of each other in the plot. We also see G3 as a separated correlated domain. 

In the second row of Fig. |l3| we give the results for realization B, in which G3 and G$ share some of their spins. 
Those spins are marked by (^)- In column (c) we see these spins inside G's- The rest of the spins of G3 are not correlated 
with them. Some of them are correlated with Gi and Gi-, and others seem to be in another domain. 

In the third row of Fig. |l3| we present the results for realization C in which G3 C Q' 3 . Here spins of G'3 seem to form 
a correlated set also over C\ , though the correlations are not high enough for it to be considered as a domain by our 
definition. 



VI. STATE OVERLAP 



We have presented a description of the system in its low T phase, relating state space behavior to the microscopic 
structure in spin space. Most of the previous literature, however, did not directly measure the microscopic features 
of the system but examined their indirect implications on other parameters, such as the widely addressed overlap 
distribution P(q). Beyond making contact with the literature, which concentrates on measuring P(q), the aim of this 
Section is two fold: (i) we show how our methods allow a useful decomposition of this function into its physically 
relevant constituent parts, and (ii) we demonstrate that our picture provides a microscopic interpretation of the 
observed P(q). To this end we focus here on Pj(q), the overlap distribution for a specific realization {J} of the bonds, 
whereas earlier works ]l^ , p0| , p4| , ff5[ presented results for the average over the disorder, P(q) = [Pj(q)]j. 

Two technical comments should be first made. First, because of overall spin-flip symmetry, the function Pj{q) is 
symmetric and hence we can limit our attention to q > 0. Second, since for most realizations |C?i| > N/2, we have 

l\l <*) 

where by Pj (q) we denote the distribution of overlaps between pairs of states |j,veC, so that we have to deal only 
with such pairs. 

A. Decomposition of Pj(q) and P(q) 

The overlap distribution for a specific realizations of the randomness, Pj(q), is expected to be the sum of two main 
parts 

Pj(q) = pKq) + Pj(<l), (33) 
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where P}(g) is the overlap distribution within a valley (and between a valley and its spin reversed counterpart), and 
Pj(q) is the overlap distribution between states that belong to two different valleys. P}(q) converges to S(\q\ — q i ea)/2 
in the thermodynamic limit, where q^A is the Edwards- Anderson order parameter, which will also be denoted as the 
"self-overlap". Pj{q) is the sum of several contributions, corresponding to different pairs of valleys. 

In the thermodynamic limit this separation is unambiguous; if two microstates [i and v are separated by a macro- 
scopic energy barrier, they belong to two different valleys and their overlap q^ v contributes to Pj(q). For finite 
systems this separation is problematic; our picture and method, however, does allow us to estimate Pj(q) or, to be 
more precise, to calculate a function Pj{q) defined below, which is a lower bound to it. In our picture, the transition 
between such pairs of microstates (that belong to two different valleys) is associated with flipping a specific set of spin 
domains. Consequently, having identified the relevant spin domains, we can identify when fj, and v belong to different 
valleys and also the level in the states' hierarchy at which they differ. 

A remaining apparent ambiguity concerns the level of the state hierarchy at which we "stop" and decide whether 
a particular pair of microstates belongs to different valleys or not. Suppose we stop the decomposition of C at some 
level n and denote by C™ the clusters obtained at this level. The overlaps obtained from pairs of microstates that 
belong to different valleys at this level are assigned to the distribution Pj' n (q), and pairs from the same valleys to 
Py n (q): 

C"C" 



P°/ n {q) = Y,PT aC h<i), (34) 

where, from Eq. (|33|) 

Pj(q) = Py\q) + P°f n (q) for q > 0. (35) 

Clearly, by going down a level further, to n + 1, some pairs that were assigned to P l j n {q) will be reassigned to 
Pj' n+1 (q), but if a pair was in Pj' n (q) it will stay in Pj' n+1 (q). This argument clearly shows that Pj ,n {q) obtained 
at any level is a lower bound to P°(q). This point is explained again below for the particular case of n = 2. 

To demonstrate how natural is the separation of eq. (|33|), we consider pairs of states \x € C\ and v G C%, i.e. pairs 
taken from state clusters that appear at the second (n = 2) level of the states' hierarchy. According to our picture 
such pairs contribute a non-vanishing part of Pj(q), which we denote by Pf lC2 (q) (= P°f 2 {q), since for n = 2 C has 
only these two sub clusters) This function, as well as its complement Pj(q) — Pf lC ' 2 (q) are presented, for T — 0.2 and 
L = 8 in Fig. |l4], for four realizations of the randomness. The figure shows clearly that the separation is natural, and 
not just an artifact of our analysis. 

For all these four realizations the spin domain Q2 is clearly identifiable and is "macroscopic" (note that this holds 



for more than 80% of the realizations, see table III). In all these cases the states fi and v belong to different valleys, 



and contribute to P°j{q)- There may be, however, pairs of states which also contribute to Pj(q), but are not included 
in Pj lC ' 2 (q). This happens when (at least) one of the state clusters €1,62 has internal structure and decomposes into 
sub-clusters (i.e. higher level valleys). Say C\ contains two such sub-clusters, C\ a ,Cxb- The overlap of a pair of states 
fi € Cia and v € Ci& contributes to Pj(q), and is not included in Pf lC2 (q); hence the latter function is a lower bound 
on the former. As discussed above in Sec. [v|, such internal structure of C\ (or C2) is associated with a spin domain G3 



structure of Pf lC2 (q)s which is discussed further below. 



(or Q' 3 ). This structre is clearly present for the realizations in Fig. 14 (a) and (d), as evident from the multi-peaked 
j 



We now generate a distribution P ClC2 (q) which is a lower bound on the contribution of P ClC2 (q) = [Pj ' (q)]j to 
the average distribution P(q). In order to assure that P ClC ' 2 (q) constitutes a lower bound to P ClC ' 2 (q), we included in 
P ClC2 (q) only contributions Pj lC2 (q) from those realizations J in which Q2 was relatively large, namely \§2\ > 0.057V. 
For the other realizations we set the contribution to the average over J to zero; hence our P ClC2 (q) is a lower bound 
to the true P ClC2 (q) (which, in turn, is a lower bound to P°(q)). In Figs. [l5| we show the distributions P(q) and 
P ClC2 (q). The data indicates that the weight in the tail for small q stays finite with increasing L (at least for this 
range of sizes), in agreement with earlier studies [^9| , p0|j2^^5|| which just measured P(q). For systems with Gaussian 
couplings P l (q) has a very small contribution at \q\ < 0.7 and P°(q) is the dominant part of P(q) in this range. For 
an Ising spin-glass with binary couplings, however, the difference between the distributions is significant and proper 
care must be taken when delicate issues, such as triviality of P(q), are investigated 
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FIG. 14. The distribution Pj(q) for four realizations of {J} at T = 0.2 in 3-D. The distribution in (a) is the same as in the 
top frame of Fig. [l^. The solid line describes P^ lC2 (g) and the dashed line plots the rest of the distribution, Pj(q) — Pj lC2 (g). 
The latter contains a large peak at q ~ 1 which is the distribution P}(q), of overlaps inside the valleys. 



B. Interpretation of Pj(q) in terms of spin domains 

Our aim is to interpret the distribution Pj(q), obtained for a particular realization, in terms of the state clusters Ci 
and spin domains Q a that were discussed in the previous sections. Before going into a detailed discussion and analysis, 
we state the interpretation that arises, for the four realizations whose Pj(q) was shown in Fig. The first of these, 
Fig. [l4| (a), corresponds to a system in which C 2 has internal structure, due to a sizeable domain Q' 3 ; its counterpart, 
C/3 is too small to have a clear signature. The size of Q'^ governs the splitting of the peak drawn with a solid line 
and also of the peak at high q (dashed line). In the systems of Fig. Il4| (b ) and (c) neither C\ nor C 2 have noticeable 
internal structure; the domains 03,0$ are microscopic. The system of |14| (d) has internal structure for both C\ and 
C 2 , induced by domains G3 and G'3, respectively. The sizes of these two domains govern the observed splitting of both 
the solid and dashed curves. 

One can associate each peak of Pj(q) with the overlaps of pairs of states that are related by flipping one or more 
of the previously identified spin domains. In this regard our interpretation resembles the RSB picture || which also 
relates the peaks of P(q) |2fJ] to overlaps between configurations in different valleys. 

To substantiate these claims and make them more precise we consider in detail the realization whose (ordered) state 
and spin data matrix is given in Fig. |l2|, and whose Pj(q) (shown in Fig. |lj (a)) is reproduced and magnified in Fig. 
Era. For this realization we clearly identified three spin domains; Q±, Q 2 and Q'^. Disregarding the splitting induced by 
G% (and C?3, if present) we identify two main peaks that dominate Pj C (q). We performed a fit of Pj C (q) to a sum of 
two Gaussians, 

PS C (q)^b 1 exp[(q~q 1 ) 2 /a 1 2 ]+b 2 e^[(q-q 2 ) 2 /a 2 2 } , (36) 

with Oi, hi and qi as fit parameters, yielding the dotted curves in the upper part of Fig. |l6|. The center of the (split) 
peak at low q is q\ and the high-g data is centered at q 2 . 

To see how these qi are related to our state clusters and spin domains, note that the overlap q^ u between states fi 
and v is related to the size of the set (defined in (pfj|)), of spins that flip when passing from state /x to v\ 
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FIG. 15. (a) The partial distribution P c ^ ( g ) for D = 3 L = 4, 5, 6, 8. It is normalized so that 2 J Q J P c ^ [q) is its weight 
in the total P(q). For clarity only a few representative error bars are shown, (b) The distribution P(q) for the same systems 
as in (a), (c) P ClC2 (q) as in (a) but for D = 4 L = 3,4, 5. (d) P(q) for the same systems as in (c). 



q»» = 1 - 2\G»v\/N . (37) 

For nearly all state pairs fi, v £ C the domain Q\ is in the state ffi; hence < 1 — so that q^ v > 2\Qi\/N — 1. 
The state pairs belong to one of two types: 

1. Pairs in which Q2 flips between f|"2 to JJ-2 or vice versa. These pairs contribute to Pj' 1 ^) = P C j lC2 {q). The 
definition of G2 yields that in most such cases G2 £ Q^u and hence 2\Qx\/N — 1 < q^ u < 1 — 2\Q 2 \/N. 

2. Pairs in which neither Q\ nor g 2 flip contribute to Pj'^g) = Pj lCl {q) + Pj 2 ° 2 (q) . For these pairs in most cases 
\Gnv\ <N- \Qi ug 2 \ and hence q^> 2(\Gx\ + \g 2 \)/N-l. 

The peak centered at q\ , is attributed to state pairs of the first type, and hence 

2\Gi\/N-l< qi <l-2\g 2 \/N (38) 

The other peak, centered at q 2 , is attributed to state pairs of the second type, and thus we expect 

q 2 > 2(10x1 + \g 2 \)/N-l (39) 

These two inequalities yield q 2 — qi > 2\g 2 \/N. Evidently, this structure of Pj{q) is completely consistent with our 
picture of spin domains that govern partition of state space into well defined clusters. By a detailed analysis [ |37| we 
have shown that the (at least) two-peaked structure of P.j{q) survives for large L. 
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FIG. 16. (Top:) The distribution Pj(q) for the same 3D realization whose data were presented in Fig. [L2l The dotted 



line is a fit to the sum of two Gaussians (see text). (Middle:) The partial distribution Pj 1 2 (q) for the same realization. The 
dotted line represents a fit to the sum of two Gaussians. (Bottom:) The difference between the two previous distributions. 



In some realizations, such as the ones that yield Figs. |lj(a) and |lj(d) Pj(q) has more peaks, since Pj 1 2 (q)(= 
Pj 2 (q)) exhibits two or more peaks; this splitting is due, as mentioned above, to spin domains G3 and G' 3 . We analyzed 



P^ lC2 (q) in the same way as we did for Pj' C (q), using the same form of fit as in Eq. (|36|). For example, in the middle 
part of Fig. |l6|, q^j and q\^ denote the centers of the two Gaussians, with G2 and G3 playing the previous roles of Gi 
and G2- 

For much larger systems, for which the state hierarchy is expected to have more than two clear levels, we expect to 
find a finer structure in P(q). It will exhibit multiple peaks, each related to different domain sizes. The heights and 
widths of the peaks are expected to be governed by the sizes of the state clusters that contribute to it which, in turn, 
are determined by the correlations between the spin domains that generating these clusters. Each of these peaks can 
be isolated and measured separately by observing the overlap of states of the corresponding clusters. 

The shape of P(q) we describe above resembles the one assumed by RSB. It is important to re-emphasize, however, 
that our P(q) was obtained for finite systems; its resemblance to the form predicted by RSB does not necessarily 
mean that the latter picture is the correct one. In fact, previous studies JlJ,|lJ,^,^6| of the link overlap (defined 
in Eq. (^)) indicate that it is trivial, which contradicts the RSB scenario, though this conclusion has been disputed 
in Refs. fl27|-p9||. In fact, our picture and results also do not appear to be consistent with RSB since we find a 
non-ultrametric state structure, as we show in Sec. VII. 



VII. ULTRAMETRICITY 



Ultrametricity is one of the main characteristics of the mean field RSB picture. Efforts to establish |16| or dismiss 
p7| the existence of ultrametricity in short range spin glasses did not yield conclusive results. We presented in Sec. IV 
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indications that w\2 , the width of the distance distribution between states from C\ and Ci , does not vanish, implying 
a non-ultrametric structure of state space. Here we look for a more direct test of ultrametricity. The main problem is 
that we can equilibrate only small systems, where ultrametricity is hindered by finite size effects. Ultrametricity is a 
statement about the geometrical properties of triangles formed by three "pure states" (or by three micro states that 
belong to different pure states). All three have to belong [E| to C, and for small systems only a small fraction of the 
realizations contain such triplets of states. 

For D = 3 at T = 0.2 we measured p, the fraction of realizations for which Q3 (or Q' 3 ) were large enough to induce two 



clearly separated peaks of Pj lfi2 (q) (see Sec. [VI BQ . We found, for L = 4, 5, 6, 8 the values p = 0.006, 0.026, 0.056, 0.090, 



respectively. At D = 4 the similar fractions, at T = 0.2 and for L = 3,4,5 are p = 0.02,0.030,0.080. Note that for 
both D = 3, 4, p increases with the size of the system. 

Our method of analysis allows us to identify the realizations that do contain such triangles of states and use 
exclusively them to investigate whether ultrametricity does or does not hold. In this way we avoid many finite size 
effects that might obscure the results. 

A set of objects with a distance measure D is ultrametric if any three objects a, (3 and 7 form an isosceles triangle, 
with the base equal to or smaller than the two equal sides. This demand can be formulated as the requirement that 
the inequality 

D aP < max{D Q7 , D Pl } . (40) 

be satisfied for all three choices of the distance placed on its left side. 

When the system is in the high T paramagnetic phase it will exhibit ultrametricity, since, as L — > 00 the probability 
distribution of distances will be P{D P „) — 8{D pu — 1/2) and all triangles will be equilateral. Similar behavior occurs 
inside a specific valley at T < T c , since for two states p and v inside the valley P(D pv ) — > 5{D pv — (1 — <7ea)/2), 
where qsA is the Edwards- Anderson order parameter. 

The non-trivial result of RSB is that the valleys themselves are ultrametric. In order to investigate this claim, we 
have to focus on triplets of states, each chosen from a different valley. For large systems with many valleys this does 
not require special care, since almost all triplets of states will belong to three different valleys. For small systems, 
however, a large fraction of the possible triplets will have at least two states from the same valley. Such triplets should 
be disregarded. 

Our way of analysis provides us with tools to examine ultrametricity for small systems. We utilize the state hierarchy 



obtained in Sec. IV to carefully choose triplets of states from different state clusters. We chose three clusters: C2, C\ a 
and C\b- The last two clusters are the "children" of C\ in the state dendrogram, i.e. C\ — C\ a UCn- According to our 
picture a triplet of states, one from each of these three clusters, belong to three different valleys, since we have to flip 
a correlated domain with a macroscopic number of spins in order to move from one cluster to another. To move from 
C2 to C\ we have to flip Q2 from configuration ^2 to configuration -ft 2. Similarly, when moving from C\ a to C\b we have 
to flip G3 from -ft~ 3 =f*|~3 "ffe) to 4J-3=4J-3 Orijfte) (see Subsection |VOj ). Due to the small sizes studied, in this paper 
we do not present any conclusive evidence that Q 3 is indeed macroscopic. However, if (in the L — > 00 limit) it is not 
macroscopic, our method predicts that there are only four valleys (determined by Q\ and Q2) and hence the the RSB 
picture clearly does not hold. 

In order to have a quantitative measure of ultrametricity we define an index K in the following manner. Let p, v 
and p be three states, so that D pv > D pp > D up . We define 



K, vp = °^ D »" . (41) 



The triangle inequality requires D vp > D pv — D pp so we have < K pvp < 1. Ultrametricity demands D pv = D pp so 
if there is ultrametricity we expect P(K) — > 5(K) as L — > 00. 

We measured P{K pvp ) for fx S C2, v G C\ a and p G Cu- We used our samples for T = 0.2; since as the temperature 
is lower and more distant from T c , the state structure should be clearer and less blurred by finite size effects. We 
measured the distribution of K for each realization, and then obtained P{K) by averaging over the disorder {J}. In 
all systems we found with high probability that K^ up = 1 exactly (see Table M).This happens when Q pv , the set of 
spins one has to flip when going from p to v 1 coincides precisely with Q pp U Q vpi the union of the two sets that are 
flipped when we go from p to p and to v. This is, however, clearly a finite size effect; as L increases the probability 
P(K = 1) decreases dramatically. Therefore we do not include this part of the distribution in our estimation of P(K). 
If this part of P(K) broadens as L increases, its exclusion cannot be achieved by simply ignoring the triangles with 
K = 1. This, however, is clearly not the case: we present in Table the probability P(0.9 < K < 1), and show that 
its increase with L is much too small to compensate for the decrease in P(K — 1). 

In order to disregard this finite size effect we truncated P(K — 1) from P(K) and renormalized to get the distribution 
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FIG. 17. The distribution P{K\K < 1) of K^ p , for fj, G C 2 , v G C\ a and p G Cit- All systems are sampled at T = 0.2. 

P ( A-|A-<1) = { ^)/^< 1 ) *<} . ,42, 

For large L we expect P(X = 1) to vanish, and P(K) will approach P(K\K < 1). The results are plotted in Fig. 
0. In Table we give the mean and variance of P(K\K < 1). Though we deal with small systems, it seems that 
P[K\K < 1) converges to a distribution with non- vanishing mean and variance, indicating breakdown of ultrametricity 
for the three valleys studied. 

Again one should address the question: do these results remain valid in the large L limit? We have to show that 
the state triplets we used, from C2, C\ a and Cu, have a finite statistical weight as L — ► 00. In Sec. we showed 
that \C-z\/N remains finite if the average correlation C12 between Gi and G2 does not approach one. From the same 
argument we conclude that if the correlation of G3 with Gi UG2 does not approach one then both |Ci a |/iV and \Cu\/N 
do not vanish and the weight of such state triplets remains finite, and the system does not exhibit ultrametricity. We 
do have evidence that the average correlation c(Gs,Gi U G2) of G3 with Gi U G2 m fact decreases as L increases, but it 
is not conclusive. 



VIII. SUMMARY AND DISCUSSION 



We have presented a new picture of the spin glass-phase in finite dimensional systems. This picture - State 
Hierarchy Induced by Correlated Spin domains (SHICS) - is consistent with numerical findings of a non-trivial overlap 
distribution []l9|,^o|,^4| and macroscopic spin domains which cost only a finite energy to flip fl^ , |l4| . Our results differ 
from the conventional interpretations |T(^|l7| of the droplet picture; nevertheless, the scenario presented in the original 
work of Fisher and Huse [^-||, and also the work of Newman and Stein |) 18 1, is of sufficient generality to allow 
consistency with our findings. 

In the spin glass phase, the system consists of macroscopic spin domains of variable sizes. Each of these domains 
flips as a coherent entity, and the flipping costs only a finite free energy. The variability in size gives rise to a 
hierarchical structure in state space. At each level in the hierarchy some state clusters split; each such splitting is 
associated with a spin-domain. The first (highest) level splitting (to C,C) is associated with the largest domain Gi', at 
the next level the two observed splittings (C — ► C\ , C% and C — ► C\ . C2 ) are related by symmetry and hence governed 
by the same, second largest domain G2- At each level, the state clusters are labeled according to the orientation of 
the corresponding domains. 

Below the second level, different spin domains are involved depending on which state cluster is being subdivided, 
e.g. G3 is the domain whose orientation splits the states in C\, while a different domain G'3 is involved in splitting C2. 
Although G3 ^ G-3, in general they may share some of their spins. The state space structure in the lower levels of 
the hierarchy has to be further investigated for larger systems. Specifically, one has to verify that G3 and G'3 do not 
vanish as L — > 00 . 

Some details of our hierarchical picture do not appear to be consistent with RSB. According to the RSB scenario, 
the states have an ultrametric structure, which implies that for any two state clusters defined at a certain level of the 
hierarchy e.g. C\ and C2, the distribution of overlaps qij between i € Ci and j € C2 should approach a delta- function 
for large L. We presented in Sec. indications that the width of the distribution P(Dij), of values in D, may not 
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vanish for L — > oo, indicating absence of ultrametricity. We also presented direct evidence for lack of ultrametricity 



in Sec. VII . However, studies on larger sizes arc needed to verify that the test which indicates lack of ultrametricity 



will still yield the sa me co nclusion as L — > oo. 

In Sections ^ and VII we demonstrated how, by separating the state space into its components, we can calculate 
various quantities using only a chosen part of this space, thus obtaining more reliable numerical results and reducing 
finite size effects. 

Clustering analysis can be applied also to other systems with a non-trivial phase space structure, i.e. which have 
several valleys which are not related by any apparent symmetry, such as random field models, see e.g. the discussion 
in Ref. ]7j , or other models with random anisotropy jl9| . It can help not only in the investigation of the macroscopic 
properties of a system, but also in understanding the micro-structure that give rise to its properties. 
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D 


L 


[\Gs\]j/N 


[c(g 3 ,Giug 2 )]j 


P(G 3 + 0) 


[\G's\]j/N 


P{Q'z + 0) 


[\G 3 ng' 3 \/\g 3 \]j 


P(g 3 + and g' 3 + 0) 


3 


4 


0.048±0.003 


0.55±0.015 


0.914(4) 


0.087±0.008 


0.834(6) 


0.23±0.019 


0.772(8) 




5 


0.046±0.003 


0.52±0.015 


0.914(4) 


0.085±0.009 


0.882(5) 


0.15±0.016 


0.818(7) 




6 


0.043±0.003 


0.48±0.015 


0.924(3) 


0.081±0.009 


0.896(4) 


0.19±0.017 


0.832(6) 




8 


0.036±0.003 


0.43±0.017 


0.905(5) 


0.076±0.010 


0.905(5) 


0.16±0.019 


0.827(8) 


4 


3 


0.045±0.003 


0.56±0.015 


0.928(3) 


0.094±0.010 


0.838(6) 


0.25±0.020 


0.782(8) 




4 


0.037±0.003 


0.48±0.015 


0.908(4) 


0.061±0.007 


0.920(3) 


0.16±0.016 


0.844(6) 




5 


0.034±0.005 


0.43±0.024 


0.84(1) 


0.072±0.014 


0.865(8) 


0.19±0.027 


0.73(1) 



TABLE IV. The size of the spin domain g 3 and g' 3 , the correlation of g 3 with C/i U <?2 and the relative part of g 3 and of 
g' 3 , which is common to both these spin domains. All results are taken for realizations where the domains concerned do not 
vanish, and we give also the probability of this to happen. All data was taken for T — 0.2. We present the average over 
these realizations { J} ± the statistical error, obtained by dividing the standard deviation by s/Ns, where N s is the number of 
realizations that contributed to each average. 



D 


L 


P(K=1) 


P(0.9 < K < 1) 


mean(_R" ) 


var(if) 


3 


4 


0.78 


0.0007 


0.385 


0.073 




5 


0.57 


0.0082 


0.426 


0.066 




6 


0.35 


0.0126 


0.447 


0.068 




8 


0.08 


0.0269 


0.476 


0.066 


4 


3 


0.74 


0.0012 


0.362 


0.068 




4 


0.38 


0.0116 


0.413 


0.067 




5 


0.10 


0.0095 


0.406 


0.061 



TABLE V. The third and fourth columns show the probability for K^ vp = 1 and 0.9 < K^p < 1, for fj, € C2, v e Ci a and 
p G Cib- The fifth and sixth columns give the mean and variance of the distribution of P(K\K < 1). All systems are sampled 
at T = 0.2. 
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